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Section  1 

Introduction  and  Summary 


The  primary  goal  of  this  project  is  to  improve  the  capability  to  identify  and  measure  surface 
waves  for  the  purpose  of  earthquake/explosion  discrimination.  We  develop  improved,  higher 
resolution  earth  and  dispersion  models.  The  models  consist  of  approximately  550  distinct  crust 
and  upper  mantle  structures,  with  surface  layering  and/or  ocean  depths  that  vary  on  a  one  degree 
grid.  There  are  a  total  of  64,800  earth  models  and  dispersion  curves,  but  the  tomographic 
inversion  is  performed  only  for  the  550  distinct  crust  and  upper  mantle  models,  with  the  shallow 
structure  constrained  by  other  information.  The  data  set  used  in  the  inversion  now  consists  of 
approximately  540,000  phase  and  group  velocity  dispersion  measurements  obtained  from  a 
variety  of  sources.  The  starting  models  for  the  inversion  are  a  modification  of  Crust  2.0  (Laske  et 
al.,  2001;  Bassin  et  al.,  2000)  over  AK135  (Kennett  et  ah,  1995).  Surface  sediments  are  defined 
using  the  global  sediment  maps  of  Laske  and  Masters  (1997),  and  ocean  bathymetry  is  defined 
using  the  Etopo5  topographic  data  set. 

Automatic  identification  of  surface  waves  at  the  International  Data  Centre  is  currently  performed 
by  narrow-band  filtering  the  data  at  several  frequencies,  and  then  comparing  the  arrival  times 
with  a  regionalized  dispersion  model.  We  have  implemented  and  tested  a  new  procedure  in 
which  we  first  phase-match  filter  the  data  and  then  apply  narrow-band  filters  to  the  compressed 
waveform  and  use  a  detection  test  similar  to  the  current  test.  This  allows  us  to  take  advantage  of 
the  improved  signal-to-noise  ratio  of  the  phase-match  filtered  waveforms,  while  retaining  the 
robustness  of  narrow-band  filtering  for  frequency-dependent  signal  identification.  After  phase- 
matched  filtering,  the  predicted  arrival  time  is  zero  at  all  frequencies,  so  we  test  to  see  if  the 
arrivals  are  within  a  time  window  similar  to  that  used  in  the  existing  test.  To  test  the  procedure, 
we  processed  the  same  data  set  using  five  degree  and  one  degree  models  with  and  without  phase- 
matched  filtering.  Detections  using  the  one  degree  model  with  phase-matched  filtering  increased 
by  40%  compared  to  the  five  degree  model  currently  in  use  at  the  IDC  without  phase-matched 
filtering. 

We  use  long-period  waveforms  from  historic  nuclear  explosions  to  assess  the  potential  of 
automated  Rayleigh  wave  travel-time  picks  to  improve  seismic  event  locations.  The  improved 
accuracy  of  locations  reported  by  Yacoub  (2000),  based  on  20-second  Rayleigh  waves,  provided 
the  impetus  for  this  work.  Although  surface  wave  arrival  times  cannot  be  measured  as  accurately 
as  body  wave  arrivals,  that  surface  waves  are  much  slower  means  that  accurate  locations  can  be 
achieved.  We  find  that  good  locations  can  be  determined  with  surface  waves,  which  could 
potentially  improve  locations  made  with  body  waves  alone  provided  that  the  paths  are  relatively 
short  and  an  accurate  dispersion  model  is  available  for  the  region.  This  may  be  especially 
important  for  small  events  with  few  total  measurements.  To  assess  the  value  of  single  surface 
wave  measurements  for  that  purpose,  we  also  estimate  the  accuracy  of  single-station  distance 
estimates. 
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Section  2 

Improved  Dispersion  Models 


2.1  Overview 

The  most  important  information  required  for  improving  surface  wave  detection  and  measurement 
is  accurate  global  dispersion  maps.  Consequently,  this  research  program  has  concentrated  on  the 
development  of  three-dimensional  velocity  models  of  the  earth’s  crust  and  upper  mantle  by 
tomographic  inversion  of  surface  wave  dispersion  measurements.  This  work  uses  an  extension  of 
the  technique  described  by  Stevens  and  McLaughlin  (2001),  in  which  a  global  earth  model  with 
149  distinct  model  types  on  a  five  degree  grid  was  developed.  That  model  is  now  being  used  for 
routine  surface  wave  identification  at  the  International  Data  Centre  (IDC).  The  technique  used  at 
the  IDC  is  to  compare  predicted  group  velocity  dispersion  curves  derived  from  these 
regionalized  models  with  measured  dispersion  curves  from  observed  surface  waves.  As 
discussed  later,  phase  velocities  derived  from  these  earth  models  can  also  be  used  to  develop 
phase-matched  filters  to  improve  signal  to  noise  ratio  and  optimize  the  detection  test. 

Development  of  improved  earth  models  and  dispersion  curves  has  proceeded  using  the  following 
approach: 

1.  The  starting  point  was  the  5  degree  IDC  149  model  set  (Stevens  and  McLaughlin,  2001) 
which  was  based  on  approximately  90,000  dispersion  measurements. 

2.  New  dispersion  measurements  were  added  to  the  data  set  and  the  same  model  set  was 
reinverted  with  the  new  data. 

3.  New  model  types  were  added  in  areas  with  increased  data  or  where  the  data  misfit 
indicated  that  new  model  types  were  required  and  a  new  inversion  was  performed  with 
the  larger  data  and  model  set. 

4.  A  procedure  was  developed  for  including  shallow  structure,  particularly  sediment 
thicknesses  and  ocean  depths,  on  a  one  degree  grid  while  inverting  for  a  limited  set  of 
distinct  models,  most  of  which  remained  on  a  five  degree  grid,  in  the  crust  and  upper 
mantle. 

5.  The  locations  and  boundaries  of  the  underlying  model  types  were  redefined,  starting  with 
the  Crust  2.0  model,  so  that  they  follow  plate  boundaries  and  other  geologic  constraints. 
The  inversion  technique  was  also  modified  to  allow  small  variations  in  Moho  depth  and 
layer  thickness  to  be  defined  on  a  one  degree  grid  for  each  model  type. 

With  our  current  approach,  the  inversion  is  performed  for  shear  velocity  structures  in 
approximately  550  distinct  crust  and  upper  mantle  model  types,  but  the  shallow  structure, 
bathymetry,  as  well  as  small  changes  in  layer  thickness  and  Moho  depth,  can  vary  on  a  one 
degree  grid  (64,800  distinct  models).  Below  the  maximum  inversion  depth  (~300  km)  the  Earth 
model  is  fixed  to  match  AK  135.  The  top  few  km  of  the  model  (consisting  of  water,  ice  and/or 
sediments)  is  fixed  and  matches  data  from  1  degree  bathymetry  maps  and  1  degree  sediment 
maps.  One  advantage  of  inverting  for  a  fixed  set  of  crust  and  upper  mantle  structures  is  that  it 
reduces  the  problem  to  a  more  manageable  size.  An  earth  model  consisting  of  16200  2  by  2 
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degree  cells  and  with  12  layers  per  cell  would  have  200,000  free  parameters.  Our  current  model 
consists  of  556  model  types,  which  reduces  the  total  number  of  free  parameters  to  8400.  Other 
advantages  of  this  type  of  inversion  are  that  it  uses  all  frequencies  (and  both  phase  and  group 
velocity)  at  once,  with  different  frequencies  resolving  different  length  scales;  the  technique  keeps 
similar  structures  consistent  with  each  other,  preventing  the  fluctuations  in  nearby  frequency 
points  common  in  more  traditional  group  velocity  inversions;  and  it  allows  Moho  depths, 
subcrustal  layer  thicknesses,  and  lateral  boundaries  between  geologic  regions  to  be  added  to  the 
model  as  fixed  a  priori  information.  The  principal  disadvantages  are  that  distinct  model  types 
must  be  chosen  carefully  and  added  sparingly  if  needed,  and  that  there  is  no  smoothness 
constraint  between  adjacent  structures.  The  tradeoff  here  is  that  without  an  explicit  continuity 
condition  we  can  incorporate  real  discontinuities  such  as  ocean/continent  boundaries  and 
boundaries  between  tectonic  regions,  but  the  inversions  will  in  some  cases  produce  adjacent 
models  with  anomalous  lateral  discontinuities. 

2.2  Data 

The  data  used  to  derive  our  earth  models  consist  of  more  than  540,000  measurements  of  surface 
wave  group  and  phase  velocity  dispersion  from  earthquakes  and  explosions.  These 
measurements  were  obtained  from  various  sources  and  are  for  frequencies  ranging  between 
0.005  and  0.1667  Hz  with  the  great  majority  of  observations  between  0.01  and  0.1  Hz.  The  data 
set  has  been  derived  from  a  variety  of  regional  and  global  studies  including  the  following:  global 
surface  wave  group  velocities  from  earthquakes  derived  using  PIDC  GSETT3  data  (Stevens  and 
McLaughlin,  1996),  augmented  with  more  recent  measurements  derived  from  PIDC  data;  surface 
wave  phase  and  group  velocity  dispersion  curves  from  underground  nuclear  test  sites  (Stevens, 
1986;  Stevens  and  McLaughlin,  1988),  calculated  from  earth  models  for  270  paths  (test  site  - 
station  combinations)  at  10  frequencies  between  0.015  and  0.06  Hz;  phase  and  group  velocity 
measurements  for  western  Asia  and  Saudi  Arabia  from  Mitchell  et  al.(  1 996)  for  12  paths  at  17 
frequencies  between  0.012  and  0.14  Hz;  the  global  phase  velocity  model  of  Ekstrom  et  al.  (1996) 
for  9  periods  between  35  and  150  seconds  calculated  for  each  5  degree  grid  block  from  a 
spherical  harmonic  expansion  of  order  1  =  40;  group  velocity  measurements  for  Eurasia  from 
Ritzwoller  et  al.(  1 996)  and  Levshin  et  al.(l  996)  for  20  frequencies  between  0.004  and  0.1  Hz 
with  500  to  5000  paths  per  frequency;  Antarctic  and  South  American  group  velocity 
measurements  from  the  University  of  Colorado  (Vdovin  et  al.,  1999;  Ritzwoller  et  al.,  1999); 
high  frequency  Eurasian  dispersion  measurements  from  University  of  Colorado  (Levshin  and 
Ritzwoller,  pers.  comm.,  2001),  and  a  large  set  of  dispersion  measurements  from  Saudi  Arabia 
provided  by  Herrmann  and  Mokhtar  at  St.  Louis  University. 

2.3  Tomographic  Inversion 

Certain  types  of  structures  can  significantly  affect  surface  wave  dispersion,  but  either  cannot  be 
resolved  adequately  by  the  data,  or  require  a  grid  finer  than  five  degrees  for  resolution.  It  is  best 
to  constrain  these  structural  features  using  independent  data.  The  effects  come  from  variations  in 
sediments  and  in  water  depths,  and  the  approximation,  by  fixed  size  cells,  of  boundaries  between 
zones,  for  example  plate  boundaries  and  coastlines.  Important  information  is  lost  when  sediment 
thicknesses  and  properties,  and  depth  of  water  columns  are  averaged  over  all  of  the  cells 
comprising  a  particular  model  type.  We  therefore  developed  a  method  for  incorporating 
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sediments,  bathymetry,  coastlines  and  zone  boundaries  determined  on  a  1  degree  grid  into  our 
inversions  while  retaining  a  relatively  small  number  of  distinct  crust  and  upper  mantle  models. 

Each  equation  in  the  tomographic  system  is  of  the  form  shown  below.  Sg  is  the  slowness  of  the 
gth  1  degree  cell;  nrij  is  the  shear  wave  velocity  of  the  jth  layer  where  j  spans  all  the  adjustable 
layers  of  all  model  types.  Xjg  is  the  length  of  ith  raypath  traveled  in  the  gth  cell,  ATj  is  the  travel 
time  residual  for  the  ith  raypath,  and  Xj  is  the  total  length  of  the  ith  raypath. 


1  dSx  ^  AT 

—  >  Am,  >  — —X,„  = — - 

x,j  ,**j  *  x, 


(i) 


Some  of  the  inversion  data  set  consists  of  model  interpretations,  such  as  the  Harvard  phase 
velocity  models.  For  these  cases  we  have  a  dispersion  value  at  each  point  rather  than  a  path 
averaged  value,  so  equation  1  becomes: 


Z  Amj 


dm) 


=  AS,. 


(2) 


where  ASg  is  the  slowness  residual. 

The  tomographic  inversion  has  the  form  shown  in  equations  1  and  2.  These  are  combined  into  a 
single  matrix  equation  AAm  =  S  where  there  is  one  column  of  matrix  A  for  each  model  layer, 
and  one  row  for  each  data  point.  Predicted  dispersion  and  partial  derivatives  of  the  model  are 
calculated  from  64,800  one  degree  cells  which  carry  the  information  about  sediments  and 
bathymetry  as  well  as  crust  and  upper  mantle  structure.  The  free  parameters,  however,  indexed 
with  j,  span  a  coarser  scale,  that  of  the  model  types.  These  equations,  plus  equations  for  damping 
and  smoothing,  are  solved  using  the  LSQR  method  (Nolet,  1987).  The  damping  and  smoothing 
equations  each  have  dimension  equal  to  the  number  of  model  layers  and  are  appended  to  the 
equations  above. 

Regularization  is  controlled  by  two  parameters,  one  for  a  vertical  smoothing  condition  applied  to 
layers  in  each  model  type,  and  a  damping  parameter.  These  have  been  applied  in  two  different 
forms.  In  the  earlier  5  degree  inversions,  a  smoothing  constraint  was  applied  that  minimized  the 
difference  in  smoothness  between  the  starting  model  and  the  final  model,  and  the  damping 
constraint  minimized  the  difference  between  the  shear  velocity  values  of  the  starting  and  final 
models  for  each  iteration.  More  recent  inversions  have  adopted  a  different  approach.  The 
inversions  now  apply  a  smoothness  condition  to  the  model,  but  allow  discontinuities  at  the  Moho 
and  at  the  sediment/crustal  or  ocean/crustal  boundaries.  Damping  now  uses  the  initial  starting 
model  derived  as  described  above  as  a  constraint  on  the  inversion,  rather  than  the  starting  model 
that  changes  between  iterations.  The  effect  of  this  is  to  stabilize  the  inversion  while  finding  a 
solution  that  matches  the  data  as  well  as  possible,  consistent  with  geophysical  constraints  and  a 
model  that  is  vertically  smooth  between  known  discontinuities. 
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2.4  Inversion  on  a  Five  Degree  Grid 

We  briefly  discuss  the  five  degree  models  that  were  developed  in  the  earlier  phases  of  this  project. 
The  five  degree  models  consist  of  2592  5x5  degree  cells,  each  one  associated  with  a  particular 
model  type.  Each  model  type  consists  of  plane  layers,  each  with  uniform  P  and  S  wave  velocities, 
density  and  Q,  with  layers  extending  to  a  depth  of  about  220  km.  As  with  the  one  degree  models, 
the  S-wave  velocities  are  treated  as  free  parameters  which  are  estimated  by  tomographic  inversion 
of  observations  of  group  and  phase  velocity  dispersion.  The  surface  wave  dispersion  observations 
are  a  large  subset  of  those  used  for  the  one  degree  inversions. 

Stevens  and  McLaughlin  (2001)  describe  development  of  a  global  earth  model  with  149  distinct 
models  on  a  five  degree  grid  that  is  now  used  for  surface  wave  identification  at  the  P1DC  and  1DC. 
The  development  used  global  tomographic  inversion  of  about  90,000  dispersion  measurements, 
using  the  Crust  5.1  model  (Mooney,  et  al.,  1998)  as  a  starting  point.  This  process  was  continued  by 
Stevens  and  Adams  (1999,  2000)  who  increased  the  number  of  data  points  to  more  than  248,000 
and  increased  the  number  of  earth  model  types  to  399. 

2.5  Methods  for  Refining  Inversions 

The  number  of  free  parameters  in  the  global  model  can  be  increased  as  the  number  of  data  points  is 
increased.  In  this  section,  we  describe  procedures  for  finding  new  model  types,  and  give  an 
example  showing  how  these  procedures  were  applied  for  one  of  the  five  degree  inversions. 

2.5.1  Finding  new  model  types. 

Addition  of  new  data  permits  us  to  increase  the  number  of  model  types  and  therefore  the  total 
number  of  free  parameters.  To  select  new  model  types  we  use  the  results  of  2D  group  velocity 
tomography  to  identify  cells  with  the  largest  apparent  differences  between  observed  and  calculated 
dispersion  (Stevens  and  Adams,  1999).  “2D  tomography”  refers  to  the  more  traditional  type  of 
tomographic  inversion  in  which  we  invert  observed  group  velocity  dispersion  measurements  over  a 
narrow  frequency  band  to  find  the  group  velocities  for  each  cell.  For  each  of  9  frequency  bands  we 
perform  a  tomographic  inversion  using  the  group  velocity  residuals  from  the  3D  tomographic 
inversion  (inversion  for  earth  structure  using  both  phase  and  group  velocity  at  all  frequencies)  as 
data,  and  perturbations  in  the  group  velocities  for  each  of  the  2592  cells  as  parameters.  Cells  with 
large  perturbations  for  several  frequency  bands  are  selected  as  possible  candidates  for  association 
with  new  model  types.  This  is  done  because  spurious  fluctuations  may  occur  between  frequencies 
in  the  2D  inversions,  but  a  consistent  trend  between  frequencies  is  a  good  indicator  of  a  real 
physical  effect.  In  addition,  spurious  perturbations  can  occur  in  areas  of  poor  resolution,  usually 
caused  by  the  effect  of  a  small  number  of  inaccurate  measurements.  Therefore,  as  an  additional 
guide,  we  compute  resolution  matrices  of  the  2D  tomographic  systems  of  equations.  Well-resolved 
cells  with  large  group  velocity  perturbations  are  given  new  model  types  and  poorly  resolved  cells 
with  large  perturbations  are  not  unless  there  is  other  supporting  evidence.  We  use  the  Lanczos 
method  of  approximating  SVDs  for  sparse  matrices  described  and  applied  to  tomography  problems 
by  Vasco  et  al.  (1999),  and  code  written  by  Berry  (1992). 
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Figure  1  shows  the  diagonal  elements  of  an  approximation  of  the  resolution  matrix  for  the  2D 
tomographic  system  using  group  velocity  measurements  for  frequencies  above  and  including  0.025 
Hz  and  below  0.0330  Hz.  This  approximation  is  made  by  computing  only  the  573  most  significant 
singular  values  of  the  SVD  out  of  a  possible  2592  (the  number  of  cells)  and  using  the  573 
corresponding  singular  vectors  to  approximate  R,  i.e.  R  =  VPVPT.  Vp  is  the  rectangular  2952x573 
submatrix  of  V  where  V  is  defined  by  the  Lanczos  decomposition  of  the  full  matrix  A  =  USVT. 


0.0  0.2  0.4  0.6  0.8  1.0 


Figure  1.  Map  of  the  diagonal  elements  of  the  resolution  matrix  for  group  velocity  measurements  in  the  frequency 
range  0.025  to  .0330  Hz  using  an  approximation  of  the  SVD  using  the  Lanczos  method  and  573  out  of  a 
possible  2592  singular  values. 


2.5.2  Selection  of  new  model  types  from  resolution  calculations. 

We  computed  the  resolution  matrices  for  the  2D  group  velocity  tomographies  for  eight  frequency 
bands:  below  0.01  Hz,  0.01  to  0.0140  Hz,  0.0140  to  0.0167  Hz,  0.0167  to  0.025  Hz,  0.025  to  0.033 
Hz,  0.033  to  0.04  Hz,  0.04  to  0.05  Hz,  0.05  to  0.067  Hz.  Those  cells  associated  with  diagonal 
components  of  resolution  matrices  greater  than  0.6  for  8  frequencies  were  given  new  model  types. 
The  0.6  value  for  a  cell  means  that  there  is  some  smearing  among  other  cells  but  that  the  true 
perturbation  of  the  cell  contributes  the  most  (60%)  to  the  tomographic  value.  The  new  types  had 
the  same  number  of  layers  and  layer  thicknesses  as  the  originating  types,  and  the  initial  parameters 
were  set  the  same.  In  this  example,  the  addition  of  these  new  model  types  increased  the  number  of 
types  from  230  to  388  and  reduced  the  weighted  standard  deviation  in  data  by  about  5%. 
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2.6  Five  Degree  Models 


Figures  2  through  5  show  the  dispersion  curves  for  our  best  five  degree  model. 


Figure  2.  50-second  phase  velocity. 


Figure  3.  50-second  group  velocity. 
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Figured  20-second  phase  velocity. 


Figure  5.  20-second  group  velocity. 
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2.7  Inversion  on  a  One  Degree  Grid 

As  discussed  previously,  for  the  one  degree  model,  the  inversion  is  performed  for  shear  velocity 
structures  in  approximately  550  distinct  crust  and  upper  mantle  model  types,  while  the  shallow 
structure,  bathymetry,  as  well  as  small  changes  in  layer  thickness  and  Moho  depth,  can  vary  on  a 
one  degree  grid  (64,800  distinct  models).  P  wave  velocities  are  constrained  via  a  constant 
Poisson’s  ratio  of  0.27,  and  density  via  Birch’s  Law  of  p  =  0.65/3  +  400  in  where  p  is  density 
and  /?is  shear  velocity  in  MKS  units.  The  layers  vary  in  thickness  (averaging  about  10  km)  and 
reach  down  to  a  depth  of  about  300  km.  Below  300  km  the  Earth  model  is  fixed  to  match 
AK135.  A  continuity  condition  is  applied  between  the  deepest  layer  used  in  the  inversion  and  the 
structure  below.  The  top  few  km  of  the  model  (consisting  of  water,  ice  and/or  sediments)  are 
fixed  and  match  data  from  1  degree  bathymetry  maps  made  by  averaging  Etopo5  5  minute 
measurements  of  topography,  and  Laske  and  Masters  (1997)  1  degree  maps  of  sediments. 

Our  current  preferred  one  degree  model  consists  of  556  model  types.  There  are  two  main 
varieties  of  model  types,  those  originating  from  the  Crust  2.0  2x2  degree  crustal  types  (Bassin  et 
al.,  2000  and  Laske  et  al.  2001)  and  those  based  on  ocean  ages  (Stevens  and  Adams,  2000).  The 
Crust  2.0  models  are  an  improvement  of  the  Mooney  et  ah  (1998)  Crust  5.1  model  revised  and 
refined  to  a  2-degree  scale.  Crust  2.0  consists  of  a  2-degree  partition  of  the  Earth  with  each  cell 
having  one  of  341  crustal  types.  Each  type  consists  of  up  to  8  layers:  ice,  water,  up  to  two 
sedimentary  layers,  3  crustal  layers  and  an  upper  mantle  layer.  Vp,  Vs,  thickness,  density  and  Q 
are  specified  for  each  layer.  The  341  crustal  types  fall  into  24  broad  categories.  Each  of  these 
categories  consists  of  several  types  with  the  same  crustal  velocities,  densities  and  Q,  but  differing 
in  thicknesses  of  the  layers,  sedimentary  structure,  ice  thickness  and  bathymetry. 

Because  many  of  the  Crust  2.0  models  differ  only  in  shallow  structure  rather  than  crustal 
velocity,  we  use  a  somewhat  different  set  of  models.  We  combine  similar  model  types  that  vary 
in  shallow  structure,  since  we  treat  those  shallow  variations  as  independent  constraints.  This  has 
the  effect  of  reducing  the  number  of  model  types  relative  to  Crust  2.0.  However,  we  also  separate 
models  in  widely  separated  regions  so  that  results  from  a  model  in  North  America,  for  example, 
do  not  affect  a  model  of  the  same  type  in  Eurasia.  This  increases  the  number  of  model  types.  We 
also  reparameterized  the  structure  of  the  oceans  based  on  ocean  age.  We  formed  17  ocean  types 
based  on  the  Muller  et  al.  (1 997)  isochron  map.  There  are  five  types  for  each  of  the  three  major 
oceans:  Pacific,  Indian  and  Atlantic.  The  five  types  match  the  age  ranges  0  to  10.9  Ma,  10.9  to 
20.1  Ma,  20.1  to  40.1  Ma,  40.1  to  83.5  Ma  and  83.5  to  180.0  Ma.  The  remaining  two  types  are 
for  the  Arctic  Sea  for  ages  0  to  40. 1  Ma,  and  40. 1  to  1 80.0  Ma.  Sedimentary  structures, 
bathymetry  and  lateral  boundaries  between  different  oceanic  types  are  included  on  a  1  degree 
level  of  detail.  The  mantle  structure  for  our  new  parameterization  comes  from  the  1997  version 
of  AK135  (Kennett  et  al .,  1995).  This  parameterization  serves  as  the  initial  model  for  a  series  of 
3D  tomographic  inversions.  Table  1  summarizes  the  new  parameterization. 
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Table  1.  Summary  of  the  new  Earth  parameterization. 


Model 

Element 

Mantle 

structure 

Vp  Vs  and 
density 

Moho 

depth 

Thicknesses 
of  3  crustal 
layers 

Crustal  Vs,Vp, 
density, 
and  Q 

Sediments 

Bathymetry 

Ice 

Vertical 

Resolution 

-20  km 
up  to  max. 
depth 
-  300km 

-1  km 

-1  km 

1-10  km 

1 00  meters 

1 00  meters 

0.5  km 

Horizontal 

Resolution 

Variable:- 
average  -  10 
degrees. 

2  degrees 

2  degrees 

Variable. 

Same  as  for 
mantle. 

1  degree 

1  degree 

2  degrees 

Parameter 

Origin 

New 

parameter¬ 

ization 

Crust  2.0 

Crust  2.0 

New  parameter¬ 
ization 

1  degree  sed.  map 
(  Laske  and 

Masters,  1997) 

EtopoS 

Crust  2.0 

Inversion 

Role 

Vs  is  free. 

Vp  is 

constrained  via 
fixed  Poisson’s 
ratio. 

Density  via 
Birch’s  law. 

Q  fixed. 

fixed 

fixed 

Vs  is  free. 

Vp  is 

constrained  via 
fixed  Poisson’s 
ratio. 

Density  via 
Birch’s  law. 

Q  fixed. 

fixed 

fixed 

fixed 

2.8  One  Degree  Models 

Figures  6  through  9  show  phase  and  group  velocity  maps  at  20  and  50  seconds  for  the  final  one 
degree  model. 


Figure  6.  50-second  period  phase  velocity  map. 
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Figure  7.  50-second  period  group  velocity  map. 


Figure  8.  20-second  period  phase  velocity  map. 
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Figure  9.  20-second  period  group  velocity  map. 


The  shear  velocity  depth  profiles  of  the  three  main  oceanic  types  after  tomographic  inversion  are 
shown  in  Figure  1 0.  Note  that  the  low  velocity  zone  at  ~80  km  is  deeper  for  younger  structures, 
as  we  would  expect,  and  that  the  structures  found  for  each  ocean  are  quite  similar  at  the 
corresponding  age  even  though  inversions  were  performed  independently. 


depth 

km 


Pacific  Ocean  Indian  Ocean  Atlantic  Ocean 


Vs  km/s 


Figure  10.  Shear  wave  velocity  depth  profiles  for  the  main  15  ocean  types.  These  results  are  consistent  with 
thermal  models  of  oceanic  lithosphere. 
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Section  3 

New  Detection  Algorithm  Using  Phase-Matched  Filtering 


Automatic  identification  of  surface  waves  at  the  International  Data  Center  is  currently  performed 
using  the  processing  program  Maxsurf  by  narrow-band  filtering  the  data  at  several  frequencies 
and  then  comparing  the  arrival  times  with  a  regionalized  dispersion  model.  An  automatic  surface 
wave  processing  program,  Maxpmf,  has  been  developed  which  is  similar  to  Maxsurf  except  that 
it  applies  phase-matched  filtering  to  seismograms  and  calculates  path-corrected  spectral 
magnitudes  in  addition  to  Ms.  Maxpmf  integrates  a  regionalized  phase  velocity  model  to  generate 
a  phase-matched  filter  which  is  used  to  compress  the  surface  wave  waveforms.  Figure  1 1  shows 
an  example  of  waveform  compression  after  application  of  phase-matched  filters  to  a  data  set. 
Note  that  the  signal  is  both  compressed  and  enhanced  relative  to  the  background  noise.  We  want 
to  take  advantage  of  the  improved  signal  to  noise  ratio  to  improve  surface  wave  detection,  and 
have  experimented  with  a  number  of  ways  to  do  this.  The  most  obvious  method,  applying  an 
STA/LTA  detector  to  the  compressed  waveform  does  not  work  as  well  as  the  current  IDC 
procedure. 

A  procedure  that  does  work  well,  however,  is  to  use  essentially  the  same  procedure  currently 
used  at  the  IDC,  but  to  apply  it  to  the  phase-matched  filtered  waveform  instead  of  the  original 
waveform.  The  allowable  time  window  has  the  same  width  as  in  the  current  procedure,  but  the 
predicted  arrival  time  is  now  zero  in  all  cases.  Figure  12  shows  the  recommended  way  to 
perform  surface  wave  identification  using  phase-matched  filtering,  which  is  to  apply  narrow- 
band  filters  to  the  phase-match  filtered  waveforms  and  the  look  for  arrival  times  near  zero. 
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Figure  11.  Data  in  the  5  km/sec  to  2  km/sec  time  window  from  an  mb  3.9  South  Pacific  earthquake  after 
transformation  to  a  KS36000  instrument  (left)  and  after  compression  by  phase-matched  filtering  (right). 
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Figure  12.  The  narrow-band  filter  detection  test  (top)  compares  measured  group  velocity  with  a  regionalized  group  velocity  model.  To  improve  detection,  we 

phase-match  filter  the  data  first,  and  then  apply  narrow-band  filters  (bottom).  The  detection  test  is  applied  for  the  same  time  interval,  but  the  time  is  now 
centered  around  zero.  This  example  is  for  an  mb  4.2  South  American  earthquake  observed  at  BDFB,  BOS  A,  and  ULM. 
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Section  4 

Tests  of  Dispersion  Models 


4.1  IMS  Data 

As  a  test  of  both  the  improvement  in  models  and  the  ability  of  phase-match  filtering  to  improve 
detection,  we  ran  Maxpmf  on  one  full  day  of  data  and  applied  both  the  current  narrow-band 
filtered  test  and  the  phase-matched  filter  test  described  above  to  this  data  set.  The  results  are 
shown  in  Figure  13.  The  five  bars  show,  respectively,  the  number  of  detections  determined  using 
IDC  5  degree  model  with  the  current  IDC  procedure,  the  number  of  detections  using  the  best  5 
degree  model,  the  number  of  detections  using  one  of  the  early  one  degree  models,  the  number  of 
detections  using  this  same  one  degree  model  with  phase-matched  filtering,  and  detections  using  a 
more  recent  one  degree  model.  Each  improved  model,  and  the  implementation  of  phase-matched 
filtering,  shows  an  improvement.  The  number  of  detections  increased  by  40%  over  current  IDC 
detections  with  phase-matched  filtering  using  the  recent  one  degree  model.  The  first  one  degree 
model  used  in  this  test  was  a  modification  of  the  earlier  5  degree  model,  and  did  not  include  the 
modifications  discussed  earlier  to  incorporate  the  Crust  2.0  improvements. 
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Figure  13.  Improvement  in  the  number  of  detections  for  a  day  of  data  with  improvements  in  the  model,  and  with 
phase-matched  filtering. 
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4.2  Nuclear  Explosion  Data 

As  a  second  test,  we  examine  the  ability  of  the  model  to  predict  dispersion  of  surface  waves  from 
explosions  at  historic  test  sites.  Part  of  the  dispersion  data  set  used  in  our  surface  wave  study 
consists  of  data  from  these  test  sites.  Specifically,  we  are  using  surface  wave  dispersion  from 
path  corrections  generated  by  analysis  of  multiple  surface  waves  from  test  sites  (Stevens,  1986) 
as  well  as  dispersion  data  measured  by  Yacoub  (personal  communication).  The  ability  to  predict 
this  data  is  a  good  consistency  check  on  the  global  earth  model.  Tables  2  and  3  show  the  average 
and  standard  deviations  of  the  group  velocity  residuals  for  5  major  test  sites,  and  shows  how  they 
have  changed  with  expansion  of  the  data  sets  and  improvement  in  the  models. 


Table  2.  40-second  group  velocity  %  average  residuals  (Standard  Deviations). 


Source 

IDC  5° 

149  Models 

5° 

388  Models 

1° 

556  Models 

NTS  (59) 

0.14(1.40) 

0.20(1.17) 

0.19(1.17) 

East  Kazakh  (40) 

0.75  (2.14) 

0.34(1.15) 

0.41  (1.00) 

Mururoa  (13) 

-0.81  (1.47) 

-0.86(1.88) 

-0.47(1.77) 

Novaya  Zemlya  (94) 

-0.31  (1.68) 

-0.37(1.53) 

0.28(1.62) 

Amchitka  (55) 

0.23  (1.36) 

0.43(1.13) 

0.44(1.18) 

Table  3.  20-second  group  velocity  %  average  residuals  (Standard  Deviations). 


Source 

IDC  5° 

149  Models 

5° 

388  Models 

1° 

556  Models 

NTS  (58) 

-0.52  (2.29) 

-0.22(1.80) 

-0.46  (2.29) 

East  Kazakh  (40) 

-0.79  (2.00) 

-0.51  (1.46) 

0.08(1.30) 

Mururoa  (11) 

0.39(1.44) 

0.50(1.57) 

-0.34(1.50) 

Novaya  Zemlya  (99) 

-0.39  (3.55) 

-0.18(2.92) 

0.01  (3.31) 

Amchitka  (54) 

0.70  (3.51) 

0.66  (3.10) 

0.22  (3.58) 

Although  both  the  best  5  degree  model  and  the  one  degree  model  are  significant  improvements 
over  the  IDC  model,  results  are  mixed  for  the  best  5  degree  model  versus  the  one  degree  model. 
This  is  surprising  considering  the  significant  improvement  in  detection  capability  discussed  in 
the  last  section  and  overall  improvement  in  data  fit.  A  number  of  the  long  oceanic  paths  are 
actually  worse  in  the  1  degree  model,  in  spite  of  the  improvement  in  accuracy  of  ocean 


16 


SAIC-01/1084 


modeling.  Consequently,  we  have  reviewed  the  surface  wave  data  to  identify  the  reasons  for  the 
remaining  data  misfit.  Specifically,  we  are  looking  for  the  following: 

1 .  Errors  in  the  data.  Some  of  the  data  was  well  constrained  by  multiple  signals  recorded  on 
digital  instruments,  however  some  of  the  measurements  are  based  on  a  single  or  small 
number  of  hand-digitized  waveforms.  Some  of  these  measurements  therefore  may  not  be 
accurate. 

2.  Consistent  differences  between  the  predicted  and  observed  dispersion  indicating  errors  in  the 
model.  This  indicates  that  the  inversion  needs  improvement,  which  can  be  accomplished  by 
improvement  of  starting  models,  or  creation  of  additional  model  types. 

3.  Paths  that  cannot  be  modeled  adequately  by  the  great  circle  approximation.  There  appear  to 
be  only  a  few  of  these,  and  they  consist  primarily  of  long,  grazing  paths  along 
continental/ocean  boundaries.  Figure  14  shows  an  example  of  waveforms  on  the  path  from 
NTS  to  MAJO.  There  are  three  distinct  arrivals  that  consistently  occur  for  each  event.  The 
second  arrival  corresponds  well  to  the  predicted  great  circle  arrival  time,  however  there  is  a 
larger  first  arrival  which  apparently  takes  a  faster  oceanic  path.  Figure  1 5  shows  the  great 
circle  path  from  NTS  to  MAJO. 


t — 1 — ' — 1 — 1 — i — 1 — 1 — 1 — 1 — i — 1 — 1 — 1 — 1 — r 


iii... _ i _ i _ i _ i _ . _ i _ i _ . _ i _ i _ i _ i _ i _ i _ . _ i _ i — i — i — iJ 

0  1000  2000  3000  4000 

Time  (seconds) 

Figure  14.  Surface  waves  recorded  at  MAJO  from  5  NTS  explosions  bandpass  filtered  near  20  seconds. 
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NTS-MAJO 


Figure  15.  The  great  circle  path  from  NTS  to  MAJO  is  a  grazing  path  along  the  northern  boundary  of  the  Pacific 
Ocean. 


In  order  to  improve  the  predictive  ability  of  the  global  models  for  these  paths,  we  recommend  the 
following  changes  for  future  work: 

1 .  Remove  or  correct  erroneous  data.  There  are  four  common  types  of  errors  that  occur. 
First,  the  measurement  may  be  inaccurate  because  of  a  poor  quality  waveform.  In  this 
case  the  measurement  cannot  be  corrected  and  should  be  removed.  Second,  the  data  may 
be  good,  but  the  measurement  inaccurate  for  any  of  a  variety  of  reasons.  In  this  case  the 
measurement  should  be  corrected  if  possible,  removed  otherwise.  Third,  in  the  case  of 
phase-velocity  the  wrong  phase-velocity  branch  may  have  been  selected.  Particularly  for 
more  distant  stations  the  phase  velocity  branches  (calculated  from  phases  that  differ  by 
multiples  of  27t)  can  be  difficult  to  determine.  The  improved  predictive  capability  of  the 
models  makes  it  possible  to  verify  or  correct  the  phase  velocity  in  some  cases.  Fourth,  the 
data  may  be  accurate  over  a  fixed  frequency  band  but  inaccurate  in  some  range.  There  are 
two  principal  cases  where  this  occurs.  First,  particularly  with  hand  digitized  data,  the 
group  velocity  cannot  be  measured  accurately,  most  commonly  at  the  lowest  or  highest 
frequencies.  Second,  in  the  case  of  oceanic  paths,  an  averaged  path  may  have  extended 
into  the  frequency  band  where  the  dispersion  curve  drops  rapidly.  This  typically  happens 
only  at  the  highest  frequencies  (>  0.05  Hz). 


18 


SAJC-01/1084 


2.  Improve  starting  models.  As  discussed  previously  we  are  using  Crust  2.0  as  a  starting 
model  and  constraint  on  the  inversion.  That  is,  the  inversion  tries  to  minimize  the  data 
misfit  together  with  the  misfit  between  the  model  and  the  constraining  model.  So  to  the 
extent  that  the  starting  models  can  be  made  more  realistic  the  inversions  will  improve. 
One  way  to  accomplish  this  is  to  use  the  inversion  results  to  replace  starting  models  in 
regions  that  are  well  constrained  by  the  data  while  retaining  Crust  2.0  in  other  regions,  or 
to  perform  more  detailed  manual  inversions  of  subsets  of  data  corresponding  to  particular 
regions  of  interest. 

3.  Remove  data  that  cannot  be  predicted  adequately  using  great  circle  paths.  As  discussed 
above,  some  surface  wave  arrivals  cannot  be  adequately  predicted  by  using  the 
approximation  that  the  surface  wave  follows  a  great  circle  path.  For  most  paths  this  is  a 
second  order  effect  because  although  the  surface  wave  follows  whatever  path  takes  the 
shortest  time,  any  faster  non-great  circle  path  will  also  be  longer.  So,  for  example,  if  the 
surface  wave  follows  a  path  that  is  5  %  faster  and  4%  longer  than  the  great  circle  path, 
the  travel  time  change  is  only  1%.  For  cases  such  as  grazing  paths  on  ocean  continent 
boundaries,  however,  the  effect  can  be  large.  It  may  be  possible  to  correct  for  these 
anomalous  paths  using  non-great  circle  ray  tracing,  however  the  predictive  capability  of 
these  algorithms  needs  to  be  evaluated,  otherwise  predictions  overall  could  be  degraded. 
The  recommended  approach  for  now  is  to  document  these  anomalous  paths  and  remove 
them  from  the  data  set. 

Some  of  the  errors  discussed  here  had  little  effect  on  earlier  inversions  because  the  results  were 
not  accurate  enough  to  resolve  these  differences,  but  now  that  dispersion  predictions  from  the 
inversion  models  are  becoming  accurate  to  within  about  one  percent  for  many  paths,  errors  of  a 
few  percent  in  the  data  can  significantly  degrade  the  results.  Further  improvement  will  therefore 
require  reviewing  the  data  set,  removing  erroneous  values,  and  adding  new  higher  quality  data 
points  where  they  are  needed. 
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Section  5 

Use  of  Surface  Waves  to  Improve  Location 


Surface  wave  arrival  times  can  be  used  in  the  same  manner  as  (or  together  with)  P-wave  arrival 
times  to  determine  source  location.  With  a  good  regionalized  group  velocity  model  and  well 
determined  group  arrival  times,  the  location  could  be  determined  quite  accurately,  and  in  cases 
where  P-wave  coverage  is  poor,  or  azimuthal  coverage  is  limited,  using  surface  wave  arrivals  has 
the  potential  to  significantly  improve  location  accuracy.  Yacoub  (2000)  found  that  he  was  able  to 
determine  the  location  of  7  NTS  explosions  more  accurately  with  surface  waves  narrow-band 
filtered  in  the  17  -  23  second  period  band  than  with  P  waves.  This  was  done  using  a  constant 
group  velocity  of  3.0  km/sec  to  determine  location  for  all  paths.  This  result  is  attributed  to  the  fact 
that  surface  wave  group  velocities  are  much  slower  than  P-wave  velocities,  so  that  errors  in 
arrival  time  correspond  to  significantly  smaller  errors  in  distance  than  the  corresponding  errors  in 
P  wave  arrival  times,  and  to  the  consistency  in  surface  wave  arrival  times. 

We  have  performed  some  experiments  using  our  global  group  velocity  models  and  dispersion 
measurements  to  assess  the  ability  to  locate  events  using  surface  wave  arrivals.  Table  4  lists  location 
errors  from  Yacoub  (2000)  derived  using  both  surface  waves  and  P  waves,  and  the  location  errors 
we  derived  from  independent  measurements  of  data  at  several  frequencies  from  the  same  events 
using  the  regionalized  dispersion  curves  discussed  earlier.  The  results  confirm  that  location 
estimates  are  somewhat  better  than  can  be  obtained  from  P-wave  arrivals  at  the  same  stations, 
although  the  location  errors  are  slightly  larger  than  obtained  by  Yacoub.  This  may  be  because 
Yacoub  used  a  larger  and  higher  quality  data  set.  Also,  it  should  be  recognized  that  Yacoub’s  P- 
wave  measurements  used  only  data  from  stations  that  also  had  surface  wave  data,  and  considerably 
better  results  could  be  obtained  from  a  larger  P-wave  data  set.  ISC  locations  for  these  events,  based 
on  all  available  P-wave  measurements,  have  errors  on  the  order  of  just  a  few  kilometers. 

Yacoub  (2000)  found  remarkably  constant  group  velocity  measurements,  all  very  close  to  3.0 
km/sec.  We  cannot  confirm  this  and  assume  that  it  must  be  due  either  to  a  particular  choice  of  paths, 
or  some  anomaly  in  the  way  the  measurements  were  made.  Figure  16  shows  the  distribution  of 
measured  group  velocities  from  all  of  the  explosion  data  from  all  test  sites.  While  the  mean  of  the 
distribution  is  very  close  to  3.0  km/sec,  there  is  considerable  variation  from  this  value  in  the  data  set. 

Table  4.  Location  errors  for  3  methods  and  the  number  of  stations  used  for  each  (  Yacoub  P  and  surface  wave 
locations  used  identical  sets  of  stations). 


Event 

#  STA 

Location  Error  (km) 

Yacoub 

IDC 

Yac.  LR 

Yac.  P 

IDC  LR 

Hearts 

19 

17 

11.2 

18.3 

22.5 

Backbeach 

19 

2 

9.7 

33.1 

NA 

Scantling 

18 

11 

25.1 

49.9 

18.9 

Lowball 

17 

12 

16.3 

28.3 

18.1 

Mizzen 

14 

9 

9.8 

32.0 

34.2 

Strake 

14 

0 

13.0 

26.0 

NA 

Sheepshcad 

12 

16 

9.2 

16.1 

13.0 

Averages 

13  ±  6 

29  ±  1 1 

21  ±7 
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Group  Velocity 

Figure  16.  Histogram  of  measured  group  velocities  of  20  second  Rayleigh  waves  from  explosions  at  all  major  test 
sites. 


50  S  Period  40  S  Period 


Figure  17.  Mean  and  two  standard  deviation  error  bounds  of  distance  estimates  from  single  stations,  for  different 
frequencies.  There  are  approximately  100  data  per  distance  bin. 
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Figure  1 7  shows  the  estimated  error  as  a  function  of  distance  for  single  station  measurements.  In 
general,  the  location  error  increases  with  distance  to  the  observing  station.  Lower  frequencies, 
however,  degrade  much  more  slowly.  However,  at  shorter  distances  higher  frequency  data  gives 
more  accurate  locations. 

Figure  1 8  shows  the  location  error  derived  from  surface  wave  measurements  from  Semipalatinsk 
test  site  explosions  as  a  function  of  the  number  of  recording  stations.  The  results  are  similar  at 
NTS,  with  a  typical  location  error  of  about  30  km.  It  is  important  to  note,  however,  that  location 
accuracy  depends  as  much  or  more  on  which  stations  recorded  the  event  as  on  how  many. 
Referring  to  Figure  1 7,  it  is  clear  that  surface  waves  recorded  at  distances  of  3000  km  or  less 
provide  much  better  constraints  on  location  than  more  distant  stations.  Surface  waves  are  most 
likely  to  be  useful  in  constraining  location  therefore  when  they  are  measured  at  shorter  distances 
and  combined  with  a  small  number  of  P-wave  measurements.  The  P-wave  measurements  on  their 
own  may  not  be  very  accurate,  but  even  a  single  surface  wave  accurately  measured  at  a  distance 
of  <2000  km  could  significantly  improve  the  location. 


Semipalatinsk  Explosions 


Figure  18.  Location  errors  for  Semipalatinsk  explosions.  The  median  error  for  events  recorded  at  6  or  more  stations 
is  30  km. 
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Section  6 

Definition  of  Ms  and  Path  Corrected  Spectral  Magnitudes 


One  of  the  principal  motivations  for  this  project  is  to  improve  the  effectiveness  of  the  Ms:mb 
discriminant  (Figure  19.  See  also  Stevens  and  Day,  1985;  Douglas  et  al,  1971;  Marshall  and 
Basham,  1 972)  by  reducing  the  threshold  for  which  surface  waves  can  be  detected  and  accurately 
measured. 


Ms:mb 


Eq  Ms 
Exp  Ms 
Ms=1 ,4'mb-3.25 


mb 


Figure  19.  Ms:mb  plot  for  a  data  set  of  earthquakes  recorded  at  the  PIDC  and  historical  explosions. 


In  addition  to  detection,  it  is  also  important  to  be  able  to  measure  the  surface  wave  accurately, 
and  it  is  important  for  the  surface  wave  magnitude  to  be  determined  in  a  way  that  is  as  free  as 
possible  from  regional  biases  and  variations  due  to  the  frequency  and  distance  at  which  it  is 
measured.  The  definition  of  Ms  adopted  by  IASPEI1  is 

Ms  =  log(  A/T)  + 1 .661og  A  +  0.3  (3) 

where  A  is  the  instrument  corrected  zero  to  peak  amplitude  in  nanometers,  A  is  the  source  to 
receiver  distance  in  degrees,  and  T  is  the  period  at  which  A  is  measured,  with  T  measured  near 


1  International  Association  for  Seismology  and  Physics  of  the  Earth’s  Interior. 
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20  seconds.  NEIC2  uses  the  IASPEI  formula  measuring  A  from  the  vertical  component  of  a  long 
period  seismogram  and  selecting  the  largest  amplitude  for  periods  between  1 8  and  22  seconds. 
The  factor  of  1 .66  in  Equation  3  is  intended  to  correct  for  a  worldwide  average  surface  wave 
attenuation  rate,  although  a  number  of  studies  (e.g.  Marshall  and  Basham,  1972)  have  shown  that 
this  in  incorrect  at  distances  less  than  30°.  Rezapour  and  Pearce  (1998)  recommended  redefining 
Ms  as 


Ms  =  log^/7’  +  AlogA  +  ylog(sinA)  +  xAloge  +  £> 


(4) 


which  is  based  on  the  theoretical  formula  for  surface  wave  attenuation,  with  k  =  1/3  and  y  =  0.0105. 
The  constant  D  =  2.484  makes  this  definition  of  Ms  equal  to  IASPEI  Ms  at  83°.  Equation  4  is  now 
used  as  the  standard  magnitude  definition  at  the  IDC. 

Ms  as  defined  in  Equations  3  or  4  can  be  quite  variable  at  different  stations  because  of 
differences  in  dispersion  and  attenuation  along  the  travel  path.  Also,  it  is  difficult  to  reliably 
measure  a  20  second  amplitude  from  a  time  domain  waveform  at  distances  less  than  about  20 
degrees.  Because  of  this,  a  number  of  authors  have  suggested  modifications  to  make  Ms  more 
consistent.  Marshall  and  Basham  (1972),  for  example,  derived  a  distance  correction  with  a 
smaller  slope  at  shorter  distances,  and  empirical  frequency  dependent  corrections  to  the 
amplitude  depending  on  the  type  of  earth  structure  between  the  source  and  receiver.  Note  that 
Equation  4  can  be  regionalized  by  varying  y  and  D  as  a  function  of  location,  and  it  can  be  made 
to  work  at  periods  other  than  20  seconds  by  making  y  and  D  frequency  dependent.  However,  any 
time  domain  magnitude  suffers  from  two  problems:  the  difficulty  of  measuring  a  time  domain 
waveform  at  a  specific  frequency,  particularly  at  regional  distances;  and  the  variability  of  the  k 
factor  in  Equation  4  which  depends  on  whether  the  measured  surface  wave  propagates  with 
normal  dispersion  or  as  an  Airy  phase. 

6.1  Spectral  Magnitudes  and  Scalar  Moment 

The  key  to  developing  a  surface  wave  magnitude  that  is  regionalizeable,  and  gives  consistent 
values  at  regional  and  teleseismic  distances  and  in  different  frequency  bands,  is  to  use  the 
equation  for  surface  waves  in  a  plane-layered  structure  to  correct  the  spectrum.  This  equation  can 
be  factored  into  functions  that  depend  on  the  source  and  receiver  earth  structure  and  the  phase 
velocity  and  attenuation  integrated  over  the  path.  The  displacement  spectrum  for  a  Rayleigh 
wave  at  distance  r  from  an  explosion  is  given  by: 


U(a>,r)  =  M0 


Sf  (o},hx)S2(a))exp[-y p((D)r  +  i(<p0  -  cor  /  cp((o))] 
Jae  sin  (r/oe) 


(5) 


depends  on  the  source  region  elastic  structure  and  the  explosion  source  depth.  S2  depends  on 
the  receiver  region  elastic  structure,  yp  is  the  attenuation  coefficient  that  depends  on  the 
attenuation  integrated  over  the  path  between  the  source  and  receiver.  cp  is  the  phase  velocity 
integrated  over  the  source  to  receiver  path,  (po  is  the  initial  phase  of  the  source.  ae  is  the  radius  of 


2  National  Earthquake  Information  Center  of  the  U.  S.  Geological  Survey. 
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*  3  p2 

the  earth.  M0  =  —  M0  where  Mo  is  the  explosion  isotropic  moment.  This  definition  is 

introduced  so  that  the  function  S *  does  not  depend  explicitly  on  the  material  properties  at  the 
source  depth. 

We  can  use  Equation  5  to  define  a  spectral  magnitude  corrected  for  distance  and  spectral  shape. 
We  define,  for  any  event,  earthquake  or  explosion,  the  path  corrected  spectral  magnitude,  or 
scalar  moment: 


M0  = 


U(a>,r  ,6)1 


S* (<*> A )S2 exp[-/  {co)r  +  i(p0  -carl c  (<y))] 


Jc/e  sin(r/ae) 


(6) 


Although  the  imaginary  part  of  the  exponential  is  removed  by  the  absolute  value,  it  is  shown 
here  explicitly  because  in  practice  the  phase  is  used  to  generate  a  phase-matched  filter  to 
compress  the  signal  and  improve  signal/noise  ratio  prior  to  taking  the  spectrum.  The  spectrum  is 
then  averaged  over  a  frequency  band  to  smooth  the  spectrum  and  get  a  stable  measurement. 

For  an  isotropic  explosion  source  at  depth  hx ,  M0  is  constant.  Equation  6  therefore  corrects 

completely  for  all  frequency  dependent  and  distance  dependent  factors.  Using  M0or  log{  M0)  as 
the  definition  of  a  new  type  of  magnitude,  the  magnitude  value  will,  ideally,  be  identical  when 
measured  at  any  distance  range  and  over  any  frequency  band.  For  an  earthquake  with  double 
couple  moment  Mq  ,  M0  is  given  by: 


M0  =  Mq\s?(<othq,ff)/S?(co,hx)\ 


(7) 


where  S?  depends  on  the  double  couple  orientation  and  depth,  takeoff  azimuth  #and  source 
region  elastic  structure.  In  general  M0  for  an  earthquake  is  not  completely  frequency 
independent,  but  it  is  partially  corrected  for  frequency  dependence  by  removal  of  the  path 
attenuation  and  receiver  structure  and  similarities  in  the  explosion  and  earthquake  excitation 
function  in  the  same  source  region.  The  remaining  differences  mean  that  the  earthquake 
magnitude  will  vary  somewhat  when  measured  over  different  frequency  bands  while  the 
explosion  will  not.  In  particular  the  spectra  of  deeper  earthquakes  will  decline  more  rapidly  with 
increasing  frequency  (This  is  a  potential  discriminant  that  is  investigated  in  the  following 
section).  By  defining  the  scalar  moment  with  Equation  6,  we  obtain  a  measure  of  surface  wave 
magnitude  that  is  independent  of  range,  nearly  independent  of  frequency,  and  regionalizeable. 

The  functions  S*  and  S2  depend  only  on  the  source  and  receiver  points  and  can  be  stored  in  a 
simple  lookup  table.  The  functions  yp  and  cp  depend  on  the  source  to  receiver  path  and  can  be 
found  by  integrating  along  a  great  circle  path  between  the  source  and  receiver  in  a  regionalized 

earth  model.  Note  that  log{  M0 )  is  equivalent  to  Ms,  except  that  the  distance  correction  is 
replaced  by  regionalized  corrections  and  the  spectral  amplitude  is  obtained  directly  from  the 
spectrum  instead  of  measuring  a  time  domain  amplitude  and  estimating  the  frequency.  Ms  and 
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log(  M0 )  are  related  approximately  by  log(  M0  )  «  Ms  +  11 .75.  Figure  20  shows  log{  M0 )  plotted 
vs.  mb  for  a  set  of  1DC  data  and  historical  explosion  data  equivalent  to  Figure  19  for  Ms:mb.  Here 
log{  M0 )  was  calculated  using  a  frequency  band  of  0.02-0.05  Hz  for  all  data. 


log(M0):mb 


mb 

Figure  20.  Log  M0:mb  plot  for  a  data  set  of  earthquakes  recorded  at  the  PIDC  and  historical  explosions. 


6.2  Regionalization 

Equation  6  requires  regionalized  source  and  receiver  functions  and  attenuation  coefficients.  While 
the  phase  velocity  is  not  required  to  calculate  scalar  moment,  it  can  also  be  regionalized  and  used 
to  construct  a  phase-matched  filter.  Complete  regionalization  therefore  consists  of  the  functions 
Sf ,  S2 ,  y  ,  c,  and  u  (group  velocity)  evaluated  and  stored  for  each  grid  block  over  the  range  of 
frequencies  of  interest.  These  functions  have  been  calculated  for  all  of  the  earth  models  developed 
under  this  project  and  are  part  of  the  data  distribution  which  is  described  in  Appendix  A. 
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The  idea  of  using  spectral  shape  as  a  discriminant  was  first  suggested  by  Tsai  and  Aki  (1971) 
based  on  observations  of  small  earthquakes  and  explosions  on  and  near  the  Nevada  Test  Site. 
The  theoretical  basis  for  this  is  Equation  6,  which  is  flat  for  explosions  but  not  flat  for 
earthquakes.  In  particular,  the  spectra  of  surface  waves  from  earthquakes  are  predicted  to  fall  off 
more  rapidly  at  higher  frequency  because  of  the  typical  earthquake’s  greater  depth.  This  effect 
also  trades  off  with  the  effectiveness  of  the  Ms:mb  discriminant  for  higher  frequency  regional 
surface  waves.  A  shallow  earthquake  will  most  likely  be  identified  as  an  earthquake  by  the 
Ms:mb  discriminant,  while  a  deeper  earthquake  may  fail  the  Ms:mb  discriminant  when  using  a 
high  frequency  Ms,  but  may  still  be  identified  as  an  earthquake  on  the  basis  of  its  spectral  shape. 
The  excitation  function  ratio.  Equation  6,  is  shown  in  Figure  21  for  several  depths  in  a  Eurasian 
earth  structure.  This  example  was  calculated  for  an  earthquake  with  strike  0,  dip  80,  and  rake  15, 
observed  at  an  azimuth  of  45  degrees.  The  explosion  is  at  the  reference  depth  of  1  km. 


Scalar  Moment  Estimate 


Frequency  (Hz) 


Figure  21.  Scalar  moment  estimate  for  an  explosion  and  earthquake  (Equation  2.5)  calculated  for  several 

earthquake  depths.  The  spectral  amplitude  decays  more  rapidly  with  increasing  frequency  for  deeper 
sources. 


The  surface  wave  spectral  shape/depth  discriminant  has  not  been  widely  used  since  it  was  first 
suggested  because  of  the  variability  of  spectral  shape,  and  dips  and  peaks  that  result  from  path 
effects  and  other  causes  not  related  to  the  source,  particularly  on  long,  teleseismic  paths.  Now 
that  we  have  improved  regionalized  earth  models,  however,  we  can  better  assess  the  viability  of 
spectral  shape  as  a  depth  discriminant. 
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Figure  22  shows  the  calculated  path  corrected  spectral  magnitudes  for  all  events  with  known 
depth  for  day  2000333  (Julian  day  333  in  year  2000).  Here  we  have  averaged  the  spectra  over  a 
narrow  range  near  each  frequency  point  to  smooth  out  spectral  dips,  averaged  the  spectra  for  all 
stations  observing  each  event,  and  divided  by  the  value  at  0.01  Hz.  The  results  resemble  the 
calculated  values  of  Figure  21  in  some  respects.  The  deepest  event  at  133  km  does  have  the 
steepest  spectral  slope,  and  the  slope  generally  increases  from  shallow  to  deep.  However  the 
depth  dependent  spectral  decrease  with  frequency  does  not  appear  to  be  as  strong  as  predicted  in 
the  data.  Based  on  these  results,  the  surface  wave  spectral  shape  does  not  appear  to  be  consistent 
enough  to  make  a  reliable  depth  discriminant. 


Scalar  Moment  Measurement 
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Figure  22.  Observed  path  corrected  spectral  magnitudes  for  events  with  a  range  of  depths  for  day  2000333. 
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Section  8 

Conclusions  and  Recommendations 


This  study  has  focused  on  development  of  improved  methods  for  detection  and  measurement  of 
surface  waves.  As  part  of  this  effort,  we  have  developed  a  gradually  improving  set  of  global 
earth  models  suitable  for  calculation  of  phase  and  group  velocity  dispersion  curves.  Over  the 
course  of  this  study,  the  earth  models  have  evolved  from  a  five  degree  model  based  on  90,000 
dispersion  measurements  to  a  one  degree  model  based  on  over  500,000  measurements.  The 
inversion  procedure  allows  shallow  structure  and  Moho  depth  to  vary  on  a  one  degree  grid  over 
approximately  550  distinct  crust  and  upper  mantle  models.  In  section  4,  we  discussed  detailed 
recommendations  for  improving  the  results  further.  We  have  used  phase-velocities  derived  from 
these  models  to  form  phase-matched  filters,  and  developed  an  improvement  on  the  surface  wave 
detection  algorithm  that  uses  phase-matched  filtering  followed  by  narrow-band  filtering  for 
detection.  We  performed  a  test  on  one  day  of  data  at  the  P1DC  with  several  dispersion  models 
and  found  a  significant  improvement  using  this  technique.  We  have  also  tested  the  use  of  surface 
waves  for  location.  Although  we  cannot  confirm  improvement  as  significant  as  suggested  by 
Yacoub  (2000),  it  does  appear  that  surface  waves  can  add  an  additional  constraint  that  could 
significantly  improve  location,  particularly  in  cases  where  there  are  only  a  few  P  arrivals  and  the 
surface  waves  are  measured  on  short  paths. 
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Appendix  A 
Data  Deliverable 


This  report  is  accompanied  by  a  data  deliverable  that  can  be  obtained  on  request  as  directed  by 
the  contracting  agency.  The  data  deliverable  is  in  the  form  of  a  Unix  tar  file  and  contains  the 
following: 

1 .  A  directory  named  “str”  containing  64800  earth  models.  The  name  of  each  file  is  constructed 
from  the  longitude,  colatitude,  and  base  structure  name. 

2.  A  directory  named  “LP”  containing  regionalized  dispersion  curves  and  other  quantities 
derived  from  the  earth  models.  Files  are  as  follows: 

LP  grid.LR  -  grid  that  identifies  one  earth  structure  with  each  one  degree  cell. 

LP_veI.LR  -  group  velocity. 

LP_pvel.LR  -  phase  velocity. 

LP_S1.LR  -  surface  wave  source  region  amplitude  function. 

LP  S2.LR  -  surface  wave  receiver  region  amplitude  function. 

LP  ellip.LR  -  ellipticity. 

LPgamma.LR  -  attenuation  coefficients. 

The  directory  also  contains  the  same  files  with  ".o"  extension  which  contain  the  grid  and  data 
file  in  binary  form  to  improve  performance. 

3.  A  file  named  “strnumbers”  listing  the  structure  name  and  the  corresponding  structure 
number  used  in  LP  grid.LR. 

4.  A  directory  named  “bin”  containing  programs  LPdisp3  and  LPdisp3_frq  compiled  under  the 
Sun  Solaris  operating  system  which  can  be  used  to  calculate  the  integrated  velocity  between 
any  two  points  on  the  earth. 

Definitions  of  SI  and  S2  may  be  found  in  Stevens  and  McLaughlin  (2001). 
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